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A high efficiency Stirling Radioisotope Generator (SRG) is being developed for possible use in long duration 
space science missions. NASA’s advanced technology goals for next generation Stirling convertors include 
increasing the Carnot efficiency and percent of Carnot efficiency. To help achieve these goals, a multi-dimensional 
Computational Fluid Dynamics (CFD) code is being developed to numerically model unsteady fluid flow and heat 
transfer phenomena of the oscillating working gas inside Stirling convertors. Simulations of the Stirling convertors 
for the SRG will help characterize the thermodynamic losses resulting from fluid flow and heat transfer between the 
working gas and solid walls. The current CFD simulation represents approximated 2-dimensional convertor 
geometry. The simulation solves the Navier Stokes equations for an ideal helium gas oscillating at low speeds. The 
current simulation results are discussed. 


Nomenclature 

T e Expansion Space Gas Temperature, K 

T h Heater Gas Temperature, K 

T mix Expansion Space/Heater Mixed Gas Temperature, K 
T r Regenerator Gas Temperature, K 

Cp Heat Capacity, J/kg-K 

m(t) Time-Dependent Mass Flow Rate, kg/s 

T(t) Time-Dependent Temperature, K 

Q Energy, W 
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I. Introduction 

High-efficiency Stirling Radioisotope Generators (SRG’s) are being developed for potential NASA space 
science missions. The SRG is being developed by the Department of Energy (DOE), Lockheed Martin (LM), 
Stirling Technology Company (STC), and NASA Glenn Research Center (GRC) for possible future multi-mission 
use, including electric power for unmanned Mars rovers and deep space missions. 1 The SRG is expected to have a 
system efficiency of about 22 percent and a specific power of about 4 W e /kg. The high SRG efficiency would reduce 
the amount of required radioisotope by a factor of four or more compared to currently-used Radioisotope 
Thermoelectric Generators. In addition, goals for advanced technology development have been established for 
second- and third-generation Stirling radioisotope power systems. Third-generation goals include achieving 
8 to 10+ W e /kg specific power and 30 to 35 percent system efficiency. 

Six Technology Demonstration Convertors (TDC’s), developed by STC as part of the SRG project, are on test at 
GRC. Operating in a dynamically-balanced opposed pair configuration, the TDC’s have demonstrated efficiencies 
from 25 to 28 percent. 

The relatively high convertor frequency and small engine clearance dimensions make detailed measurements 
difficult. Dynamic gas temperature and velocity measurements are not recorded whereas dynamic pressure and 
piston and displacer motion are recorded. Data acquisition accompanied by other means of gaining knowledge of the 
internal processes, such as using CFD prediction codes, can help identify areas for performance improvements. One- 
dimensional (1-D) design codes do not rigorously model the effects of manifold regions where the working gas 
encounters abrupt changes in flow area. It is expected that an accurate multi-D code will help to better understand 
the flow and heat transfer patterns and also help quantify thermodynamic losses throughout the internal working gas 
volumes. Such codes will be used to computationally determine flow and heat transfer for a given complex 
geometry describing the internal gas/solid interfaces of the Stirling convertor. The improved understanding gained 
from the development and experimental validation of the multi-D code will also be used to further improve the 1-D 
design codes. 

A. Stirling Computational Codes 

NASA GRC has been involved in developing and/or validating 1-D Stirling codes since the late 1970’s. 2 1-D 
codes are now used for the design of Stirling convertors. HFAST, a harmonic code, was acquired by GRC from 
Mechanical Technology Inc. (MTI) to conduct engine predictions for the 12.5-kWe Component Test Power 
Convertor (CTPC). 3 GLIMPS, developed by Gedeon Associates, was also used to compare losses in the CTPC and 
later was developed by Gedeon Associates into the Sage code. 15 Both GLIMPS and Sage 15 are one-dimensional 
design and optimization codes and have been used as the primary design tool by several Stirling engine/cooler 
manufacturers. The 1-D codes tend to over predict power and/or efficiency estimates by 10 to 20 percent, because 
they over-simplify flow and heat transfer processes. Furthermore, the 1-D codes do not rigorously model the losses 
that occur due to manifolds and abrupt area changes encountered in typical Stirling convertors. For example, 
uniform flow is assumed, and flow separation is not modeled. A better understanding of how to properly account for 
these flow and heat transfer processes will enable improvements to the 1-D codes. 

NASA conducted 2-D simulations with compressibility and grid deformation capabilities in the early 90’ s using 
the Computer Aided Simulation of Turbulent Flows (CAST) code, developed by Peric and Scheuerer in 1989. The 
code was made available to CSU (Ibrahim) who is now rewriting the code, known today as CSU 2-D , to include 
unsteady boundary conditions, various turbulence models, and a revised thermal non-equilibrium porous media 
model. The code will still use the finite volume method and the SIMPLE (pressure- velocity coupling) 3 algorithm. 

An approach to overcome the limitations of collecting data from actual convertors is to build specialized test rigs 
which can isolate the effects of Stirling processes. 4 The University of Minnesota (UMN) has collected flow 
visualization, velocity, temperature, and eddy transport data from such devices under grant work for GRC. An 
important outcome of the collaborative effort between CFD work performed by CSU 17 and experimental work 
performed by UMN 18 is an analytical technique for modeling unsteady flow and heat transfer in Stirling engine 
heater and cooler tubes which has been added to the 1-D Sage code by Gedeon Associates. 

One commercially developed computational fluid dynamics (CFD) code being used currently at GRC, CFD- 
ACE+, is available from ESI Group (formally CFD Research Corporation). CFD-ACE+, used to model the current 
simulation NASA-2D , uses the SIMPLE algorithm in solving Navier-Stokes equations (or Reynolds Averaged 
Navier-Stokes equations when turbulence is assumed). A Fluent code license has recently been purchased as well. 
The purpose of the modeling effort is to accurately model flow and heat transfer phenomena present inside Stirling 
convertors. This includes quantifying and comparing loss mechanisms associated with heat transfer and flow 
through various manifolds and clearance gap regions. Such simulations must accurately model the geometry 
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describing the internal gas/solid interfaces of the Stirling convertor in order to quantify the thermodynamic losses 
present. Understanding and reducing such losses should help achieve the goal of improved system efficiency. 

II. Stirling Engine Thermodynamics 

Stirling engines are reciprocating machines that operate on a pressure wave produced by a temperature 
differential and the motion of the power piston and gas displacing piston (displacer). Helium is commonly used as a 
working gas due to high thermal conductivity, low density which reduces pressure drop, and is preferred to 
hydrogen due to the difficulties of maintaining operating pressure with hydrogen, which can leak through solid 
walls. Reciprocating machines, in general, often contain components that move in a manner adequately 
approximated by a sinusoidal function. Typical Stirling engine components are shown in Fig. 1. The oscillating flow 
is single phase, low Mach number, zero mean flow. 4 The power piston and displacer oscillate out of phase, which 
results in a time-dependent variation of volume and thus requires use of the compressible flow equations for 
simulation despite the very low Mach numbers. This variable volume causes a fluctuation in the pressure about the 
mean operating pressure. The operating frequency is typically 50 to 105 Hz and, thus, the expansion and 
compression phases last only a few milliseconds. 4,8 
The temperature fluctuations in the expansion and 
compression spaces, mainly due to inadequate time 
for equilibrium heat transfer, are the main causes of 
adiabatic losses. 8 

The out-of-phase motions of the piston and 
displacer affect the instantaneous velocities 
throughout the internal working spaces. The flow 
undergoes acceleration and deceleration at different 
spatial locations at the same time in a cycle 
depending on the direction and position of the 
moving components. UMN has conducted velocity, 
temperature, and eddy transport measurements of 
laminar/turbulent transition for oscillating flow in 
pipes. 7 More recently UMN has collected velocity 
and flow visualization measurements of oscillating 
flow through a 90 degree turn test section and 
through a large-scale regenerator matrix. Their 
results suggest that transition from laminar to 
turbulent flow and back to laminar flow may occur 
in various regions of the engine over each engine 
cycle. It is likely that the flow is laminar for one 
spatial location and turbulent for another spatial 
location at the same moment in time due to the out- 
of-phase variable volumes. The use of low- 
Reynolds turbulence models, such as the two- 
equation k-co turbulence model, may be applicable. 

The k-co turbulence model has shown good 
agreement 7 with UMN’s 90 degree turn test section flow visualization, velocity measurements, and heat transfer 
measurement data. Also, a low Reynolds k-s turbulence model has also shown good agreement with Kornhauser’s 
experimental data for compressible non-acoustic flow in a completely enclosed volume with a moving boundary. 3 

Multi-dimensional, time-dependent simulations of the unsteady conduction and convection heat transfer 
encountered during the thermodynamic cycle does increase the cost of conducting such simulations. The first 
approach for solving this problem is running such simulations on a 64-bit, 32-processor cluster, currently being 
developed at NASA GRC. In addition to such computational resources, more efficient methods are being 
investigated for improved speed and accuracy. Such Ultra High Fidelity methods, being investigated by Dyson, 5 
resolve a simulation containing increased information on a coarse grid domain, resulting in a highly accurate 
solution with up to possibly a 1000 times increase in solution acceleration. 


Stirling Engine Components 

1 Expansion Space 

2 Heater | 

3 Displacer 

4 Regenerator 

5 Appendix Gap I 

6 Seal 

7 Cooler 

8 Compression Space , 

9 Power Piston Face 


Geometry Not To Scale 
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Figure 1. — Schematic of Stirling Engine Components 
Interfacing Working Spaces. 
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III. Stirling Thermodynamic Losses 

Thermodynamic 2 nd law available energy 
losses (or entropy increases) in Stirling devices 
are of the following general types: 4 (1) mass flow 
across a pressure difference (viscous dissipation), 
(2) heat transfer across a temperature difference, 
and (3) mixing of fluids at different initial 
temperatures. As previously mentioned, 1-D 
codes do not rigorously model the losses that 
occur due to manifolds and abrupt area changes 
encountered in typical Stirling convertors. Due to 
the complexity of calculating various local and 
instantaneous losses throughout the computational 
domain, parasitic losses are typically 
independently evaluated by analytical techniques 
and added into the 1-D simulation. The loss terms 
are often approximations based on simplifications 
of the complex geometry and flow conditions. 
CSU (Ebiana) is currently leading an effort to 
incorporate 2 nd law analysis post processing into a 
CFD-ACE+ simulation of MIT test rigs that 
experimentally simulated some of the Stirling 
device internal fluid flow and heat transfer 
processes. 16 The various parasitic losses 
associated with flow and heat transfer throughout 
the working volumes captured in the current 
modeling effort (NASA-2 D) are listed along with 
other energy flow paths in Fig. 2 and discussed in 
the following sections. 


Energy Flow for NASA-2D 

1 Energy Added 

2 Heater and Reg. Viscous 
Dissipation Loss 

3 Cooler and Reg. Viscous 
Dissipation Loss 

4 Expansion Space Adiabatic Loss 

5 Expansion Space Hysteresis Loss 

6 Enthalpy Flux and Gas 
Conduction Loss 

7 Appendix Gap Shuttle Heat 
Transfer and Hysteresis Loss 

8 Appendix Gap Pumping Loss 

9 Wall Conduction 

10 Output Work 

11 Energy Rejected 



3 IU 




Geometry Not To Scale 

Figure 2. — Energy Flow for NASA-2D CFD Simulation. 


A. Heat Transfer Losses: Mixing, Reheating, and Hysteresis 


1. Working Spaces 

Temperature fluctuations in the expansion and compression spaces result from inadequate time for heat transfer 
and cause adiabatic losses. 8,9,11 The adiabatic losses can be large in regions where the heat transfer rate is low 
between the gas and solid walls (i.e., a “nearly” adiabatic process). One component of the adiabatic loss, mixing, 
occurs during the expansion stage of the cycle. As the total working volume is increased, cooling occurs due to the 
expansion of the gas resulting in a lower expansion-space gas temperature, T e , than had entered. Under these 
conditions, gas exiting the heater then enters the 

Expansion 
Compression 


expansion space at temperature T h > T e . The 
mixing of these two gases results in a mixed 
gas temperature still lower than the temperature 
at which it had entered the expansion space (T h 
> Tmix). Work potential (or available energy) is 
decreased due to this mixing of gases at 
different temperatures. The other, often larger, 
component of the adiabatic loss, reheating, 
occurs when the gas at T mix reenters the heat 
exchanger and needs more energy from the heat 
exchanger to return to the higher gas 
temperature, T h . Figure 3 shows time-dependent 
volume-weighted cell temperature integrated 
over all volumes contributing to the expansion- 
space gas and compression- space gas as 
predicted by NASA-2D’s first cycle. The 
discontinuity (variable difference at the 


NASA-2D Working Space Temperature 

Expansion Space 



Time, s 


Figure 3. — NASA-2D Expansion-Space and Compression- 
Space Gas Temperatures. 
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beginning and end of the cycle) indicates that the simulation solution needs to be performed for many more cycles to 
reach periodic steady state. The expansion and compression spaces are not the only regions subject to adiabatic 
losses. Any low surface-to-volume ratio 4 (nearly adiabatic) gas passages adjacent to any of the components that are 
nearly isothermal promote such losses. 

Losses associated with transient heat transfer, such as hysteresis losses, can be large if the heat transfer processes 
occurring in the working spaces are neither isothermal nor adiabatic. 4,10 Hysteresis losses are due to heat transfer 
across the temperature difference between the working space gas and the surrounding walls (i.e., pressure vessel and 
displacer walls). Gas springs are devices with heat transfer processes bearing some similarity to those occurring in 
the expansion and compression spaces in Stirling convertors; however, gas springs have no external flow into or out 
of the variable volume. Adiabatic gas springs experience no hysteresis loss because there is no heat transfer across 
the temperature difference between the gas and wall. Isothermal gas springs experience no hysteresis loss because 
there is no temperature difference between the gas and the wall. Hysteresis loss for a gas spring is maximized at 
some operating condition intermediate between the two idealized extremes of adiabatic and isothermal operation. 
Hysteresis losses also occur in Stirling expansion and compression spaces which experience large temperature 
fluctuations over the cycle and heat transfer between the gas and surrounding solid walls (across a temperature 
difference). Some hysteresis loss occurs in all Stirling working space components due to gas temperature 
fluctuations over the engine cycle and heat exchange with the surrounding walls. 

The reheat losses and hysteresis losses both fall into the previously mentioned category of thermodynamic losses 
due to heat transfer across a temperature difference, so these would presumably be lumped together in a component- 
by-component 2 nd law of thermodynamics analysis, such as that being developed by CSU (Ebiana). 

2. Regenerator 

The regenerator may be the most loss-sensitive component in the Stirling convertor. The regenerator is subject to 
several heat transfer related losses that affect the performance commonly referred to as regenerator efficiency. 
Hysteresis losses in the regenerator should be small due to the “almost isothermal” regenerator gas and matrix 
temperatures at each axial location. The predominate loss in the regenerator is due to “reheating” and “recooling” of 
the gas as it enters and exits the regenerator. During the portion of the cycle when gas flows from the hot end toward 
the cold end, hot gas at a temperature approximating the heater temperature (T h ) is pumped into the matrix, where it 
heats the solid material. This energy addition to the solid matrix results in a lower regenerator gas temperature (T r < 
T h ). Similarly, the gas temperature in the cooler is lower than the regenerator gas temperature 
(Tr > Tc) due to imperfect heat transfer (i.e., a finite heat transfer coefficient and non-zero temperature difference). 
During the portion of the cycle when gas flows from the cold end to the hot end, the regenerator gas is pumped into 
the heater again where, in a situation of perfect heat transfer, the gas temperature exiting the regenerator would be 
equal to the heater temperature (T r = T h ). This is not the case in practice, where the resulting gas temperature 
entering the heater is lower than the heater temperature. Energy must now be added to the gas exiting the 
regenerator to obtain approximately the heater temperature (T h ) before it moves into the expansion space. This 
reheating requirement in the heater and “re-cooling” requirement in the cooler, both due to imperfect regenerator 
heat transfer, lead to a net enthalpy flux loss through the regenerator from the hot to the cold end. Radial heat 
transfer from the regenerator matrix and gas to the surrounding walls is typically neglected in 1-D code calculations. 
Such heat transfer can be accounted for in multi-dimensional simulations. 

3. Appendix Gap 

The appendix gap, the gas domain between the displacer and displacer cylinder, is subject to a relatively 
complicated combination of losses that needs to be better understood. Shuttle loss is assumed to take place when the 
displacer wall, at top dead center position, collects energy via conduction and convection. Then, at bottom dead 
center position, the displacer wall rejects energy to the lower-temperature gas in the appendix gap, creating an 
overall higher gas temperature in that time and spatial location. Some of that energy is added to the displacer 
cylinder wall and is then conducted to the cold side. The secondary effect of having lost the energy via shuttle heat 
transfer is the loss of work potential, or hysteresis loss. There also may be appendix gap losses related to hot gas 
from the expansion space being pumped into the appendix gap and losing heat to the wall (referred to as “pumping” 
loss). Also, the appendix gap plus displacer seal form a high flow-resistance flow path in parallel with the primary 
flow path through the heat exchangers. The effect of this secondary flow path is not very well understood. 

4. Environment 

There are also heat transfer losses from the Stirling convertor to the surrounding environment. These are a 
combination of conduction, convection, and radiation. The current CFD simulation, NASA-2D , does not consider 
losses to the surrounding environment and, instead, defines all exterior boundaries to be adiabatic. 
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B. Conduction Losses 

1. Wall and Gas Conduction Paths 

Solid conduction losses are generally well predicted analytically in 1-D codes. They are not included in the 
conservative partial differential equations but are rather handled separately and added at some point in the 
simulation. The task becomes increasingly difficult for 1-D codes to handle conduction paths when a model is 
expected to account for all conduction paths for a complex geometry. CFD codes inherently consider all conduction 
paths throughout modeled components. However, conduction through the gas is more uncertain for 1-D and multi-D 
codes, particularly in the regenerator, due to the effect of thermal dispersion resulting from eddies formed by flow 
across wires. This increases the effective gas conductivity (above the gas molecular conductivity). 

The conduction paths present in the NASA-2D simulation consider conduction in the axial and radial directions 
for solids and fluids. The following lists contain the conduction heat transfer paths modeled and omitted from the 
current simulations. 

Modeled: The following conduction heat transfer paths have been included in NASA-2D. 

• Displacer dome and body walls region (wall thickness equal to one thermal penetration depth) 

• Seal and appendix gap gas passage 

• Displacer cylinder wall 

• Pressure vessel wall (adiabatic at outer surface) 

• Regenerator porous matrix solid material (thermal equilibrium between gas and solid assumed) 

• Regenerator porous matrix gas volume (thermal equilibrium between gas and solid assumed) 

• Compression-space flow passage wall or “spider” region (wall thickness equal to one thermal 
penetration depth) 

• Piston face region (wall thickness equal to one thermal penetration depth) 

Omitted: The following conduction heat transfer paths have not been included in NASA-2D. 

• Regenerator thermal dispersion (thermal equilibrium between gas and solid assumed) 

• Piston cylinder wall 

• Displacer rod wall (wall thickness equal to one thermal penetration depth) 

• Displacer rod seal gas passage 

• Piston seal gas passage 

• Displacer interior gas volume (flexures, shields, etc.) 

• Bounce Space (large variable volume on the alternator side of the piston face) 

C. Radiation Losses 

Radiation heat transfer is possible but not modeled at this time in the current CFD modeling effort. The radiation 
heat transfer path of most interest at this time is internal to the displacer dome. 

D. Viscous Dissipation Losses (mass flow across a pressure difference) 

7. Regenerator 

The largest viscous dissipation loss is usually in the regenerator, largely due to the small hydraulic diameter of 
the porous material. 

2. Working Spaces 

Flow through a sudden change in area experiences viscous interactions at solid-fluid interfaces and produces 
pressure drops, affects heat transfer, and tends to produce flow separations which produce mixing effects. 4,6 These 
area changes include interfaces between the cooler and compression space and between the heater and expansion 
space as well as transitions via a plenum between the cooler and regenerator and regenerator and heater. 1-D codes 
use empirical expressions to represent the impact of these “end effects” on pressure and heat transfer. 
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IV. NASA GRC Stirling CFD Simulation (NASA-2D) 

The current CFD simulations were 
developed based on modifications of 
simulations performed by CSU under 
Grant NAG3-2482. The current GRC 
simulation geometry is a 2-D 
approximation of an actual 3-D free- 
piston Stirling engine geometry 
representative of the TDC’s on test at 
GRC (Fig. 4). For the NASA-2D 
simulation, heater and cooler heat 
exchanger flow passages and 
compression-space flow passages have 
been approximated by concentrically 
aligned annular channels. In order to 
conserve the flow and heat transfer 
characteristics of the 3-D heat 
exchangers, the 2-D approximations 
maintained length and inner/outer heat 
exchanger diameters while closely 
matching hydraulic diameter, heat 
transfer area, and flow volume. 

Additional geometric simplifications 
include using 90 degree corners instead of the actual finite radii on the displacer body and compression space flow 
passages. The approximated 2-D geometry, shown in Fig. 5, was then revolved about the axis of symmetry one 
degree of revolution, resulting in a 3-D model. The model needed to be revolved by a thickness of at least one grid 
cell so that the “arbitrary interface” option in CFD-ACE+ could be used (available for 3-D only) to represent the 
sliding interaction between stationary and moving meshes. The simulation (. NASA-2D ) does not actually contain 
w-velocity (circumferential velocity vector) components and, therefore, is essentially a 2-D axi-symmetric 
simulation. 

The reference pressure and operating frequency used in the simulation are 27 bar and 80 Hz, respectively. Piston 
and displacer motions are specified and therefore are not affected by the developing thermodynamic solution (where 
mass, momentum, and energy are conserved). 



Figure 4. — TDC’s on Test at GRC. 
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A. Volume, Boundary, and Initial Conditions 

The CFD model represents the working gas enclosed by solid walls. The grid applied is structured Cartesian with 
a constant cell number for all volumes (no adapting grids in deforming volumes). The oscillating piston and 
displacer are represented by moving walls and result in instantaneous pressure and temperature fluctuations, 
requiring the use of varying gas properties. The properties 
for the working gas, and solid materials are used as inputs 
to the code definition. The working gas is ideal helium for 
which the Prandtl number is set equal to 0.707; this allows 
the thermal conductivity to become a function of heat 
capacity and dynamic viscosity. The specific heat of the 
working gas is held constant. The dynamic viscosity is 
allowed to vary by employing Sutherland’s Viscosity Law, 
which corresponds to using the Sutherland curve fits for 
laminar viscosity. The density is allowed to vary as a 
function of temperature and pressure, according to the 
ideal gas equation for helium. The heat capacity and 
thermal conductivity for the pressure vessel material 
(IN718) are given as polynomial functions of temperature, 
while the density is held constant. The volumes 
representing heat exchanger fins for the heater and cooler 
have been omitted. The heat exchanger fin surfaces are 
held at constant temperatures of 650 °C and 80 to 120 °C 
respectively over the cycle. To speed up the initial steady 
solution (before piston and displacer begin to move), 
volumes throughout the entire grid domain were set to 
initial temperatures. As the steady solution matured, the 
simulation temperatures eventually converged upon an 
axial temperature difference dominated by the heat 
exchanger temperatures. This simulates a Stirling engine 
with no component motion as energy is added and rejected 
(flow and heat transfer “modules” are used so gas convection can occur during this steady solution). The unsteady 
time-dependent solution then uses the steady solution as initial conditions. The initial steady solution provides 
temperature gradients along the length of the regenerator porous medium as well as the pressure vessel and displacer 
walls, and significantly reduces overall computational time. Figure 6 shows boundary conditions for adiabatic and 
symmetry (used for axisymmetric modeling) boundary conditions and the fluid/solid volume conditions are defined. 

B. Laminar/Turbulent Flow 

So far, only laminar flow has been assumed. Several Reynolds Averaged Navier-Stokes (RANS) turbulence 
models are available in CFD-ACE+. However, the simulation currently must use either laminar or turbulent 
equations throughout the entire solution domain. This contrasts with the actual flow situation, as now understood 
based on oscillating flow rig testing, where flow can transition from laminar to turbulent and back over the cycle, 7 
and it is likely the flow is laminar for one spatial location and turbulent for another spatial location at the same 
moment in time during the cycle. A better understanding of laminar/turbulent transition in time and space is needed 
to properly model such flow physics. 

C. Porous Media Modeling 

The regenerator contains a porous material that absorbs energy during one half of the cycle and releases energy 
during the other half. The porous medium accomplishes this by having a large surface-to-volume ratio and good 
thermal conductivity. Conduction in the axial direction is undesirable, while conduction along the radial direction is 
desirable. The current “porous media module” available in CFD-ACE+ is based on a thermal equilibrium solution 
that assumes the solid and fluid are at the same temperature at any point along the length in the axial direction. The 
axial temperature distribution is generally linear and largely varies due to the working temperatures of the hot and 
cold ends of the engine. Solid-gas thermal equilibrium seems to be a very poor assumption for modeling oscillating 
flow through Stirling regenerators, since it appears to rule out the possibility of accurately representing the enthalpy 
flux loss through the regenerator. One possible alternative would be to use the porous media module to represent just 
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Figure 6. — Model Boundary Conditions. 
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the gas phase and the “user scalar” module to represent the solid phase, so that solid and fluid could be at different 
temperatures. However, this approach requires the assumption of a heat transfer correlation to represent the heat 
transfer between the solid and the gas. A DOE funded research effort recently explored regenerator matrix data 
collection (UMN) and modeling (CSU). 19 

D. Temporal and Spatial Discretization 

The time step size used in the transient simulation plays a critical role in the solution convergence and accuracy. 
An appropriate time step must be employed to compute an accurate solution. The most compressed grid positions 
coincide with the points at which the displacer reaches top dead center and bottom dead center and the piston 
reaches top dead center (alternator bounce space not modeled). Near these most compressed positions, more 
iterations and a smaller time step are required for convergence. A user-defined time step routine has been developed 
to control what time step size is used as time progresses during the computation of the thermodynamic cycle. 
Significant solution acceleration (2 times) can be achieved by appropriately increasing and decreasing the time step 
size as opposed to simply running the entire solution at the small time step size required during the compressed grid 
positions of the cycle. 

An equally important component of model accuracy is the grid mesh within the solution domain, or spatial 
discretization. The mesh plays a critical role in the calculation of an appropriate resolution of physical waves that 
best represent the actual flow physics. 12,13,14 The backward Euler time discretization scheme (1 st order accurate) is 
currently being used, while the central differencing (2 nd order accurate) spatial scheme is used for spatial 
discretization. The current mesh contains almost 110,000 cells. Evolution to a 3-D engine geometry will greatly 
increase this number (possibly by a factor of 500). Simulation convergence for the steady time dependence was 
based on monitoring pressure and temperature in particular regions of the working gas volumes. These variables 
were monitored per iteration to verify convergence. 


V. NASA-2D Results 

The energy balance was performed by calculating overall heat transfer and total PV power (addition of positive 
expansion-space power and negative compression-space power) for the system. Energy flowing across solid- to -gas 
interfaces was collected in the form of heat rate and then integrated over the cycle using Simpson’s rule. The 
interfaces considered were broken up into six different surface sets representative of different areas of heat transfer 
in the engine. Heat is transferred into the system through three surface sets including the heater fin walls (constant 
temperature surfaces), expansion- space walls, and regenerator-heater interfaces. Heat is transferred out of the system 
through three surface sets including the cooler fin walls (constant temperature surfaces), compression-space walls, 
and regenerator-cooler interfaces. Equation (1) shows the net heat rate for heat in and heat out of the system. 


NetHeatRate = 



& 


Q*dt 


Out 


( 1 ) 


The total PV power was calculated by collecting pressure and volume as a function of time for the expansion-space 
and compression-space. Pressure, density, and medium mass (helium in this case) were integrated over each domain 
volume. Time dependent volume was obtained by dividing mass by density at each time step. The pressure and 
volume were then integrated over the cycle using Simpson’s rule. 


TotalPVPower = 


cfp *dv*dt 

- Up *dv*dt 

L J J 

Exp J 


Comp 


( 2 ) 


The indicated efficiency is then calculated as shown in Eq. (3). 

WorkOut TotalPVPower 


effind =■ 


Heatln 




dt 


In 


( 3 ) 
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Enthalpy flux through each heat 
exchanger flow interface (not shown) 
was calculated by collecting spatially 
integrated data at each time step 
throughout the cycle and integrating over 
the cycle. Mass flow rate ( m ) was 
integrated across fluid- to -fluid interfaces. 
Temperature (T) was collected similarly 
and both were multiplied by heat 
capacity (constant over the cycle) for 
each flow interface of interest and 
integrated over the cycle, as seen in 
Eq. (4). 

(Jm*Cp*T*dr (4) 

The simulation results being shared are 
time-dependent computational solution 
results for the first cycle. The net 
indicated power seems to be roughly the 
correct magnitude, but the discontinuity 
(variable difference at the beginning and 
end of the cycle) indicates that the 
simulation solution needs to be 
performed for many more cycles to reach 
periodic steady state. Start up affects, 
boxed in Fig. 8(a), can be seen at the 
beginning of the cycle. The results 
representing periodic steady state 
operation (mature time-dependent 
solution) will eventually be compared to 
Sage results. Figure 8 shows plots of 
expansion-space and compression-space 
volume and pressure variation as 
predicted by NASA-2D. The PV 
diagrams in Fig. 8(a) show the 
expansion-space and compression-space 
pressure as a function of volume 
variation. The expansion- space and 
compression-space power, integrated 
over the cycle, are 242.5 watts and 
-116.0 watts, respectively. The time- 
dependent volume variation for each of 
the working spaces, Fig. 8(b), is 
specified by a sine wave via definition of 
piston and displacer motion in the code. 
Figure 8(c) shows the expansion- space 
and compression-space pressure over the 
cycle. The compression-space pressure 
does deviate slightly from the sinusoidal 
trend at a time that coincides with 
displacer bottom-dead-center position. 
Such phenomena (identified in red 
circle) are being explored to determine 
their origin. Two suspect causes are 
numerical dissipation due to compressed 
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Figure 8. — NASA-2D (a) Expansion/Compression-Space 
Pressure vs. Volume, ( b ) Expansion/Compression-Space Time- 
Dependent Volume Variation, and (c) Expansion/Compression- 
Space Time Dependent Pressure Variation. 
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grids and physical phenomena of the 
compression-space gas at that displacer 
position during the cycle. 

It can be seen that the average 
operating pressure is near 29 bar, not the 
27 bar that was set as the reference 
pressure at the beginning of the first 
cycle. The initial pressure will have to 
be adjusted to achieve the desired 
operating pressure of 27 bar due to the 
initial conditions not corresponding to 
one of the two mean pressure points over 
a periodic steady state cycle. Figure 9 
shows time-dependent expansion- space 
and compression-space temperature and 
heat transfer rate variation as predicted 
by NASA-2D. The results for temperature 
are based on a cell volume temperature 
integrated over all volumes contributing 
to the expansion and compression 
spaces. For the expansion- space 
temperature and compression-space 
temperature, shown in Fig. 9(a) and Fig. 
9(b), peak-to-peak magnitudes are nearly 
60K and 25 K, respectively. It can be 
seen that the temperature at the 
beginning of the cycle does not match 
the temperature at the end of the cycle. 
This discontinuity indicates that the 
simulation solution needs to be 
performed for many more cycles to reach 
period steady state. The cell area heat 
transfer rate was integrated across each 
solid-to-fluid interface and collected at 
each time step over the cycle, seen in 
Fig. 9(c). The interfaces considered were 
broken up into six different surface sets 
representative of different areas of heat 
transfer in the engine. The solid-to-fluid 
interfaces were considered separately to 
compare how much heat transfer 
contribution each of the surface sets has 
with the working gas. It was found for 
this solution, the heat transfer 
contributed by the cooler-regenerator 
interfaces (solid wall-to-regenerator 
interfaces with no plenum) accounted for 
over 20 percent of the total heat rejected. 
Approximately the same was found for 
the percentage of heat transferred by the 
heater-regenerator interfaces for total 
heat added. Heat transfer rates (based on 
this premature time-dependent solution) 
shown in Fig. 9(c) do not agree with 
what is expected from a valid, fully- 
converged periodic steady state solution. 


NASA-2D Working Space Temperature 

Expansion Space 


Expansion 
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Compression Space 
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NASA-2D Heat Transfer Rate 

for 6 Sets of Solid-to-Fluid Interfacing Surfaces 



Time, s 

Figure 9. — NASA-2D (a) Expansion-Space Time-Dependent 
Temperature Variation, (b) Compression-Space Time-Dependent 
Temperature Variation, and (c) Time-Dependent Heat Transfer 
Rate Across Solid-to-Fluid Interfaces for Six Surface Sets. 
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As the computational solution continues to mature, the energy balance will be performed to track convergence. At 
this time, the heat addition over the cycle from the heater is 170 watts, and the heat rejection from the cooler is -268 
watts, thus indicating negative power. 


VI. Concluding Remarks 

At this point in time, the results for NASA-2D are still maturing. It is believed that the energy balance is not 
resolved because the time-dependent solution is still developing and has not reached periodic steady state. The heat 
addition does not equal the sum of the indicated power and the heat rejected, as it should. The PV power does 
appear to be the correct order of magnitude for each of the working spaces. The Sage 1-D design code results are 
being used as a benchmark to understand how well the prediction capability of the multi-D code is proceeding. It is 
anticipated that a more mature multi-D code and its associated validation work will aid in improvements to 1-D 
design codes and to Stirling convertor design. The multi-D code needs improvement in the porous media modeling, 
turbulence modeling, and time-to-convergence. These issues are being addressed through in-house computational 
work at GRC and through university grants. GRC is installing a parallel computer cluster that will greatly enhance 
the capabilities of obtaining these types of solutions in a reasonable amount of time. GRC is also investigating new 
time efficient numerical techniques to speed up simulation convergence. 5 
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